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We consider the effect of evaporation on the aggregation of a number of elastic objects due 
to a liquid’s surface tension. In particular, we consider an array of spring-block elements 
in which the gaps between blocks are filled by thin liquid films that evaporate during the 
course of an experiment. Using lubrication theory to account for the fluid flow within the 
gaps, we study the dynamics of aggregation. We find that a non-zero evaporation rate 
causes the elements to aggregate more quickly and, indeed, to contact within finite time. 
However, we also show that the number of elements within each cluster decreases as the 
evaporation rate increases. We explain these results quantitatively by comparison with 
the corresponding two-body problem and discuss their relevance for controlling pattern 
formation in elastocapillary systems. 


1. Introduction 

The aggregation of many wet hairs into a series of clumps is familiar from everyday 
examples including wet paint brushes or eyelashes wetted by tears (see Bico et a7(||2004 


Kim fc Mahadevan|2006 for example). In these scenarios, the surface tension of a liquid 


acts to minimize the liquid surface area but is resisted by the bending stiffness of the 
hairs involved. The result is elastocapillary aggregation with a number of clumps of hairs, 
separated by larger gaps. At a microscopic scale, a similar phenomenon is observed in 
the manufacture of microelectromechanical systems (MEMS) where long thin elements 
are etched (using photolithography) and then rinsed ( [Tanaka, Morigami &: Atoda||1993 ). 
In the rinse step, the elastic elements are vulnerable to the formation of clumps that can, 
if the clumps are large enough, lead to fracture (see figure [^). 

While, in many applications, elastocapillary aggregation is an undesirable feature of the 
process and is best avoided, in other situations it is exploited to form controlled patterns 


(see 


Chakrapani et al. 2004 Pokroy et al. 2009 for example). Indeed, the clumps of carbon 


nanotubes that form when a nanotube ‘forest’ is wetted and the liquid evaporated (see 
figu re [T^) have sparked considerable interest as a design strategy at microscopic scales 
de Voider fc Hart||2013 de Voider et aL|[2M3 for reviews). Despite the importance 


(see 


of other forces (including van der Waals and electrostatic forces) at these small scales, it 
is clear that capillary forces from the liquid play a vital role: aggregation only happens 
if the forest has been wetted, and varying the surface tension coefficient changes the 
properties of the pattern ( ]de Voider &: Hart||2013 Tanaka et aZ.||1993 ). 

Most previous theoretical approaches aim to understand the pattern formation in these 
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Figure 1. Elastocapillary aggregation driven by evaporation of a volatile liquid, (a) Scanning 
electron micrograph of a pattern created to be part of a MEMS device but da maged by surface 
tension forces in the rinse step. (Image reproduced from Tanaka et a/.||1993| Copyright (1993) 
The Japan Society of Applied Physics.) (b) The formation of cellular fo ams from the elastocap- 
illary collapse of a forest of carbon nanotubes. (Image reproduced from Chakrapani et al.|2004[ 
Copyright (2004) the National Academy of Sciences), (c) The simplified model considered here 
in which rigid blocks, connected to their initial position by linear springs, replace flexible beams. 
Lubrication theory in the intervening liquid gaps leads to a second order differential equation 
for the pressure within each gap, Pj{x, t), which is solved subject to a no-flux condition at a; = 0 
and imposed capillary pressure at the meniscus, x = Xj(t), (see Appendix [A|). 


systems by focussing on energy minimization arguments (see 

Bico et a/.||2004| |Py et al. 

2007 Duprat et al. 

2012 

de Voider et al. 

2011 Jung et al.\ 

2014 for example). This 


approach neglects the dynamic manner in which the pattern forms and so does not 
address any features of the pattern that might be controlled dynamically. |Boudaoud| 


et al. (2007) developed a toy ballistic aggregation model that gives predictions for the 


distribution of cluster sizes. However, only recently have detailed dynamic models been 
developed to capture the interaction between fluid flow and elasticity that occurs in this 
problem. Several early authors developed models that couple the elastic deflection with 


the fluid flow between a pair of elastic beams (Aristoff et al. 2011 Duprat et al. 2011 


Taroni & Vella 20121. In particular, Taroni & Vella (20121 showed via their detailed model 


that even the simple system consisting of a fixed volume of liquid confined between two 
elastic beams generally exhibits multiple stable equilibria. Taroni & Vella (2012) found 
that which of these equilibria is attained dynamically depends on which is ‘closer’ to 
the initial condition, rather than which configuration has the lower energy. From this 
observation, one might expect the final state exhibited by the many-body problem to 
also be sensitive to the history of deformation, including the initial condition. 

Recently, there has been a focus on understanding dynamic aggregation in many-body 
systems, though to make progress it has been necessary to simplify the elastic problem. 


Gat & Gharib (2013) and Singh, Lister fc Vella| (|2014 ) did this by considering arrays of 


spring-block elements while Wei fc Mahadevaji[ ( |2014 ) and Wei et al. (2015) considered 
rigid rods that rotate, resisted by a torsional spring. These models are sufficiently simple 
that they allow for the development of analytical understanding of the underlying prob¬ 
lem while also being able to replicate some of the more complicated behaviour observed 
in experiments and numerical simulations. For example, Singh et al. (2014) showed that 
multiple equilibria exist, as in the full problem (Taroni fc Vella||2012 1, and that the dis- 
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tribution of cluster size predicted numerically is similar to that observed experimentally 
by Boudaoud et al. (2007). Similarly, Munro (2014) demonstrated novel scaling laws for 


meniscus motion in elastocapillary rise that resemble those found experimentally (Duprat 

erz][mi| ). 


Linear stability analysis of the simplified models of dynamic aggregation have all re¬ 
vealed that pairwise aggregation is the mode of aggregation with the largest growth rate 
(Gat & Gharib 2013 Singh et al. 2014 Wei et ar||2015 ). However, simulations and ex¬ 
periments show that larger clusters form. This apparent discrepancy is believed to be 
resolved either as the result of iterated clustering (in which quadruples form from a pair 
of pairs etc., as suggested by Gat & Gharib 2013 Wei et al. 2015) or by the propaga¬ 


tion of a disturbance from a localized perturbation, which leaves behind clusters with a 
well-defined size (as suggested by Singh et'^|2014 ). 


The models of elastocapillary aggregation discussed above have common features, most 
notably that the amount of liquid is conserved or is decreased quasi-statically. However, 
in many applications the liquid is volatile and therefore only wets the elastic elements 
transiently. Once the liquid has evaporated, elements can only ‘stick’ together if they 
were brought sufficiently close together while wet that dry adhesive forces (e.g. van der 
Waals forces) come into play. A key question is then: how close to each other do the 
elastic objects come during evaporation? In this paper we study a reduced model of this 
dynamic aggregation accounting for the effect of evaporation. 


The fluid mechanics of evaporation has seen a great deal of interest in recent years. 


with particular emphasis being placed on the evaporation of droplets (Cazabat & Guena 


2010). Much of this work was motivated by the observation that the dark ring left when 
a drop of coffee evaporates, the so-called “coffee-ring effect”, is caused by the pinning 
of the contact line and a singularity in the evaporative flux close to the contact line 
( Deegan et qZ.|[T 9^ Popov 2005). In other cases, the contact line is not pinned so that 
the evolution of both the radius and contact angle of the droplet need to be determined 
dynamically ( McHale et a/.|[2005 Stauber et a/.|[2M4 ). A variety of different evaporation 
models are available depending on, among other factors, whether the rate of evaporation 
is limited by diffusion in the vapour phase ( Murisic fc Kondlc||2011 ), conduction of heat 
from a solid boundary (Dunn et al. 2009) or convection above the drop (Kelly-Zion 


et al. 2013). However, even very simple models of evaporation give rise to extremely 


intricate fluid mechanical problems whose solution can be difficult to distinguish from 
experimental observations (Oliver et a7||2015). 


The coffee-ring effect itself may be viewed as the relatively dilute limit of the drying of a 
colloidal suspension ( Routh|2013 ). However, at higher concentrations, the densely packed 
layer of particles left after evaporation cracks just as mud cracks when it dries (iLecocq & 


Vandewalle 2002 Dufresne et al. 2006). Gracks may only propagate in these systems when 


the elastic energy that is released by fracture is sufficient to pay the associated surface 


energy penalty, much like a crack in a more classical elastic medium (Dufresne et al. 


2006). This problem therefore contains similar ingredients to elastocapillary aggregation. 


namely elasticity, capillary-induced flow and evaporation. 


In this paper, we consider the elastocapillary aggregation problem accounting for the 
motion of individual elastic elements within the system, their displacement relative to 
their initial position, and the fluid flow between elements that controls the dynamics of 
the process. We shall focus on the effect that the rate of evaporation has on the pattern 
that is ultimately formed. 
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2. Governing eqnations 

To mimic the elasticity of the beams typical in many-body elastocapillary aggregation, 
we consider an array of rigid blocks, each connected to their initial position by a linear 
spring of stiffness k (see figure [^). The gap between the blocks is considered to be filled 
with liquid of initial volume Vq (per unit depth into the page); the initial thickness of 
each gap is w. 

The surface tension of the liquid, 7 , causes a capillary suction beneath each menis¬ 
cus, which acts to bring the blocks together. This aggregation is resisted by the springs, 
mimicking the balance between capillarity and elasticity common to elastocapillary ag¬ 
gregation in many systems. 

To model the dynamics of aggregation, we assume that the gap is relatively thin, 
Vojw'^ ^ 1. In this case, lubrication theory can be used to determine the pressure profile 
Pj{x) within each gap, subject to the boundary conditions that there be no flux through 
the base (x = 0 ) and that the pressure is equal to the capillary suction at the meniscus, 
X = Xj. More details of this calculation are given in Appendix]^ the key result is that 
the pressure profile within each gap, Pj{x), is determined in terms of the values of Xj, hj 
and the rate at which the gap is growing/shrinking, dhj/At. The pressure profile Pj{x) 
can then be integrated along the wetted length of each side of the gap to give the force 
on the two neighbouring blocks exerted by the liquid in the j**' gap: 


_ dhj ^ 

^ dt hj' 


( 2 . 1 ) 


Here the gap thickness and position of the meniscus have been non-dimensionalized by 
w and Vq/w respectively, force has been non-dimensionalized by 2jVo jw'^ and time by 
the capillary-viscous time 


tr/i! — 




( 2 . 2 ) 


In (2.1) the first term represents the lubrication force provided by the liquid in the 
gap (which resists the motion if dhj/dt < 0 ) while the second term represents the net 
attractive force that results from capillary suction. 

Now, each block has two forces acting on it (from the liquid gaps on either side) and 
the difference of these must be balanced by the net spring force on the block. Subtracting 
these net forces for the and (j — 1 )®* gaps, we find that 


2K{h,-l)=F,+,-2F,+F,_,. 


Here the dimensionless spring stiffness 


K = 


kw^ 

47 Vb' 


(2.3) 


(2.4) 


To close the problem we must also give an equation for the evolution of the position of 
the meniscus, Xj, within the gap. Previous models have assumed that the total amount 


presence of convection in the atmosphere, the curvature of the meniscus and evaporation- 
induced Marangoni effects among others. Furthermore, in the two-dimensional geometry 
of interest here, the classic approach of considering evaporation limited by steady-state 
diffusion is ill-posed (essentially because of the logarithmic term in the solution of the 
two-dimensional Laplace equation). For simplicity, therefore, we assume that evaporation 
occurs at a (constant) dimensional rate e per unit surface area. It follows that the volume 


of liquid within each gap is conserved, that is d{xjhj)/dt = 0 (Gat fc Gharib|20T3 Singh 


et al. polll ). The dynamics of evaporation are, in general, complicated depending on the 
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of liquid within each liquid-filled gap evolves according to 

^{Xjhj) =-Ehj, (2.5) 

where the dimensionless evaporation rate E is defined by 


E = 



( 2 . 6 ) 


Our non-dimensionalization of the problem ensures that a uniform initial condition may 
be written as a;j(0) = hj{Q) = 1 for each liquid gap, i.e. for each j. The problem therefore 
has only two dimensionless parameters: the spring stiffness, K, and the evaporation rate, 
E. Ultimately, we shall be interested in how the final pattern, particularly the maximum 
number of blocks in any cluster, depends on these two parameters for systems consisting 
of a large number of spring-block elements. However, to gain some understanding of the 
behaviour of the system in different regions of (E, K) parameter space, we first consider 
the simpler problem of a pair of spring-block elements with a single fluid-filled gap 
between them. 


3. An isolated pair: The two-body problem 

For a pair of spring-block elements, there is only a single liquid gap, _) = 1, to consider. 


In this case Fq = F 2 = 0 so that (2.31 simplifies considerably. Denoting the (only) gap 


thickness by h(t) and the meniscus position by x(f) we find that 

x'^ dh X ,, 


dx 

dt 


= -E- 


dh X 

dt h 


(3.1) 

(3.2) 


with initial conditions a;(0) = h{Q) = 1. 


3.1. Phase-plane analysis 


The system of equations (3.1|-(3.2) is autonomous so it is natural to examine the be¬ 


haviour of its phase plane. To facilitate this, we introduce a rescaled gap aspect ratio, 
^(t) = x(t)/[Kh(t)] and a time-like coordinate s such that d/ds = K'^h^^d/dt with 
s{t = 0) = 0. The system (3.1|-(3.2) can then be written 


dh 

ds 

ds 


= hii-h-0, 

= -e [KEe + 2{l-h-0] ■ 


(3.3) 

(3.4) 


From (3.3)-(3.4) we see that this (slightly opaque) transformation has eliminated indi¬ 
vidual occurrences of the dimensionless parameters K and E] instead only their product 
KE appears. The dynamics does still depend on the separate values of K and E, however, 
since the initial conditions are now ^(0) = 1/A, h(0) = 1. Nevertheless, to understand 
the phase plane, this rescaling, and the associated reduction of the number of effective 
parameters, is extremely useful. 


It is clear from (3.3)-(3.4) that the (non-trivial) nullclines in the (^,h) plane are 


h = 1—^ and h = l—f^-\-KEff /2. Combining these with the trivial nullclines ^ = 0, h = 0, 


we find that if KE < 1/2 then the dynamical system (3.3)-(3.4) has two stable fixed 
points: ^ = 0, /i = 1 (corresponding to two dry blocks in their undeformed, separated 
configuration) and ^ = [1 -|- (1 — 2KE)^^‘^]/KE, h = 0 (corresponding to the two 
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(a) 



(b) 




Figure 2. P hase planes and numerically computed trajectories showing the behaviour of the 
system (j3.3[ )-( 3.4 1 for KE — 0.49 (a and 6) and KE = 1 (c and d). In (a) and (c) the phase 
planes show various trajectories (solid curves), corresponding to changing the value of K with 
KE fixed, with arrows indicating the direction of increasing time-like coordinate s (defined such 
that d/ds = K^h^^A/dt). The nullclines h = 1 — ^ and h = 1 — ^ KE^^/2 are shown by the 
dotted straight line and dashed curve, respectively. For KE ^1/2 there are stable equilibria at 
(0,1) (separated blocks, indicated by an open circle) and (^*,0) (stuck blocks, indicated by a 
filled circle), and two saddle points at (0,0) (indicated by -I-) and (Cs,0) (indicated by x). For 
KE >1/2 there is only the equilibrium point at (0,1) (open circle) and the saddle at (0,0) 
(-I-). (b) Trajectories of x{t) (solid curves) and h{t) (dashed curves) for K = 5 (labelled A) and 
K = 0.5 (labelled B); in each case KE — 0.49 and we observe sticking only for sufficiently small 
K. (d) Trajectories of x{t) (solid curves) and h{t) (dashed curves) for K = 5 (labelled C) and 
K = 0.5 (labelled D); in each case KE = 1 and sticking is never observed. 


blocks coming into contact as the liquid between them evaporates). We refer to the fixed 
point (^*, 0) as the ‘stuck’ configuration since we expect short-ranged forces to maintain 
contact even after all of the liquid has evaporated. There are also two saddle points for 
KE < 1/2: one at (^s, 0) with = [1 — (1 — 2KE)^^'^]/KE and another at (0, 0). 

A typical phase plane for KE < 1/2 is shown in figure [^o) (with KE = 0.49); 
the evolution of the physical variables for two example trajectories are shown in figure 
The phase plane, figurea), shows that for ^(0) = 1/K sufficiently small, i.e. for 
sufficiently stiff springs, the system ultimately tends to the separated fixed point, (0,1). 
However, for ^(0) = IjK larger than some threshold, i.e. for sufficiently weak springs, the 
system instead tends to the stuck fixed point, (^», 0). This is illustrated by the evolution of 
the physical variables shown in figure [^&). The critical value of K at which the transition 
from separated to sticking occurs, Kcnt{KE), can be determined numerically from the 
stable manifold of the saddle point at (Cs,0) (see Appendix]^. For the case KE = 1/2 
we find that the system tends to the stuck fixed point if K < 3.30 (i.e. E > 0.15) but 
otherwise becomes separated. In the limit KE <C 1 we find that iTcrit = x ~ 0{E‘^^^), 
which agrees with the result of Singh et al. (2014) in the absence of evaporation. 
















Figure 3. Regime diagrams showing the behaviour of the system as functions of the spring 
stiffness K and evaporation rate E. (a) The regime diagram for the two-body problem in which 
only stuck and separated states exist. The solid curve shows the numerically determined bound¬ 
ary between stuck and separated while the results with E — 0, R'crit = 27/4 (dotted vertical 
line), and Kait = 1/2E (dashed line), are also shown. (For K < 3.30, the boundary between 
stuck and separated states is KE = 1/2, see text.) (b) Results of numerical simulations with 
N = 2000 blocks with randomized initial perturbations showing the largest number of blocks 
in a cluster, A^max, at various {K, E). Here green shows situations in which no blocks aggregate 
ultimately while red and blue are used to distinguish results at neighbouring points of [K, E) 
space. The dashed line again shows the result TFcrit = 1/2E as a hint that the two-body problem 
is relevant to understanding the many-body problem. 


For KE > 1 / 2 , is complex so that the only stable fixed point is at ^ = 0,/i = 1 ; 
there is also a single saddle point at (0,0). Physically, we deduce that, for sufficiently fast 
evaporation or sufficiently strong springs, the system must always end up with the blocks 
separated and all of the liquid evaporated — the blocks are ‘unstuck’. An example of the 
trajectories of the system, plotted in {^,h) space, is shown in figure [^c); as expected, 
this shows that all trajectories (i.e. regardless of the value of K) ultimately reach the 
unstuck state. This is confirmed by the physical evolution along two example trajectories, 
as shown in figure]^ d). Note from the case K = 0.5, E = 2 (labelled ‘D’ in figurej^d)) 
that once the blocks start to separate they may do so very rapidly: both evaporation 
and the widening gap reduce the hydraulic resistance to motion while the total capillary 
suction in the gap is also dramatically decreased by a decrease in x. 

In summary, our analysis of the phase plane shows that for KE > 1/2 the blocks 
must ultimately tend to the separated state, while for KE ^ 1/2 the blocks stick if 
K ^ KcritiKE) and otherwise separate. This behaviour is summarized in figure]^ a), 
which shows the regions of (AT, E) parameter space for which the blocks are ultimately 
stuck or separated. However, we note that in either case the system tends to a state in 
which all of the liquid has evaporated, x = 0, whether that be with h = 1 (the separated 
state) or with h = 0, xfh = (the stuck state). 


3.2. Contact occurs in finite time or not at all 
In the limit of no evaporation, A = 0, Singh et al. (20141 showed that if the springs 


are sufficiently weak for a pair to stick, then the gap thickness h decreases with time 
according to /i ~ {3t)~^/^ for t ^ 1. In this idealized system, therefore, the blocks do not 
actually contact one another (though in reality van der Waals attraction or the roughness 
of the blocks would ensure that contact does eventually occur). 

The situation with even a small amount of evaporation, A > 0, is qualitatively different 
in that contact occurs within a finite time. To see this we let ^ 1, a;//i « in (3.1), 
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yielding 


Ah 

dt 




E 


Crucially, this is constant so that the gap thickness decreases at a constant rate and, 
hence, the gap must close and the blocks touch. 

Note that in the limit of small evaporation, E 1, and with 0 < K < Kcvit{KE) so 
that sticking occurs, we have « 2/KE leading to 

dt 

We therefore expect the time at which contact occurs to diverge as £1 —>■ 0, consistent 
with the previous observation that contact does not occur in the absence of evaporation 


(3.6) 


(Singh et al. 20141. This expectation may also be confirmed for the case of no springs. 


K = 0, where a simple calculation shows that contact occurs at t = with 


^ _ 1 ^ 2 

“ (E £13/2(2-£1)1/2 


tan 


2 

E ~ ^ 


1 / 2 ' 


-^E-^/^asE- 

V2 


0 . 


(3.7) 


4. The many-body problem 

4.1. Numerical results 

Our primary interest is in the dynamics of aggregation for a large number of spring-block 
elements. To gain some insight into the behaviour of these systems, we solved the system 


of differential equations (2.3|-(2.5| with up to Nt = 2000 elements. We used symmetry 


boundary conditions at the edge of the system, i.e. E_i = Ei and Fn^+i = Fnt-i (Singh 
et al. ]|M4l ), though similar results were obtained with free edges, Fq = = 0. Using 

symmetry boundary conditions has the advantage that the uniform state hj = 1 is an 
equilibrium. The initial condition was then taken to be a small random perturbation to 
this state, i.e. 

£,(0) = l + e7^, (4.1) 

with TZj a random number drawn from /V(0,1). For the results presented here, we took 
e = 10“^. The initial meniscus position was a;j(0) = l/£j(0) so that each gap contained 
the same volume of liquid initially. 


Two adjustments to the codes used by Singh et al. (2014) were made to account for 


the ultimate disappearance of liquid inherent in evaporation. Firstly, the liquid level 
Xj(t) was prevented from becoming negative by enforcing Xj(t) ^ 10“for all times. 
Secondly, a short range van-der-Waals interaction was introduced to ensure that once 
elements come sufficiently close to contact they remain close, but do not pass through 
one another; this mimics the sticking upon contact that is observed in MEMS devices. 
Simulations were run until \^j\ < 10“^°, so that each gap has reached an equilibrium 
to within numerical error. 

Figure shows space-time diagrams for the evolution of around 100 block positions 
during aggregation. Results are shown for a range of values of the dimensionless evapora¬ 
tion rate E but with K = 0.1 fixed. Figure]^ a) (in the absence of evaporation) shows that 
clusters first form and then shrink down as time increases. Introducing a small amount of 
evaporation (figure [^), we see that the typical size of clusters is broadly similar but that 
the clusters seem to shrink to contact within a finite time. We also observe that there 
seem to be more scenarios in which clusters start to form but ultimately split — what we 
term a ‘cluster bifurcation’. (These cluster bifurcations are also associated with a ‘kink’ 
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Figure 4. Space-time diagrams showing the position of « 100 blocks within a larger simulation 
of the elastoeapillary aggregation of Nt = 500 blocks. The spring stiffness if = 0.1 in each case 
and the evaporation rate E is varied: (a) E = 0, (b) E = 0.1, (c) E = 1, (d) E = 3. 

in the space-time diagram because the force balance is radically altered as the cluster 
splits: the new sub-clusters shift to positions where all of the spring forces are reduced.) 
Increasing the evaporation rate still further (figures we see, unsurprisingly, that 

the clusters that form shrink to contact more quickly than with lower evaporation rates. 
We also see that there are many more cluster bifurcations and that the final clusters tend 
to be smaller than those observed with smaller evaporation rates. 

Another perspective on the effect of evaporation is provided by examining the statistics 
of the cluster size N as E varies with K fixed. The inset of figure shows that the mean 
cluster size (iV) is a decreasing function of the evaporation rate, as should be expected 
from the earlier observations of figure]^ Figurej^itself shows the probability distribution 
function for a range of evaporation rates, 0 ^ i? < 30, calculated by averaging the results 
of thirty runs at each evaporation rate. To allow for a comparison of the shape of the 
distribution as the evaporation rate varies, we have normalized the a;-axis in figure by 
(TV). In these normalized coordinates, the typical (normalized) width of the distribution 
does not change noticeably as the evaporation rate increases (i.e. the width of the true 
distribution scales with the mean (N) and hence does become narrower as E increases 
and (TV) decreases). Surprisingly the behaviour in the tails of the distribution seems to 
follow quite closely that found in the absence of evaporation (this is true even at A = 30 
where only pairs of blocks or single blocks are observed). However, we note that the peak 
around TV/(TV) « 0.8 appears to become more pronounced as the evaporation rate E 
increases up to A = 1; at the largest values of E the discrete nature of the distribution 
and the decrease in (TV) blunts the peak. 

In summary, the simulation results shown in figures and suggest that the presence 
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Figure 5. The statistical properties of the cluster size N for a range of evaporation rates E, 
with K = 0.0137 fixed. Inset: The mean final cluster size, {N), as estimated, for each value of 
E, from 30 numerical experiments with initial conditions hj = 1 + eTZj, where e = 10“^ and TZj 
are random numbers drawn from A'’(0,1). The vertical error bars represent the standard error. 
The mean cluster size with no evaporation {E = 0) is shown by the horizontal dashed line. 
Main figure: The probability distribution function for the cluster size as found in the numerical 
experiments. Different evaporation rates are indicated as follows: E — 0 (►), E = 10“^ (◄), 
E — 0.5 (a), E = 1 (T), E = 5 (■) and E = 30 (•). Horizontal error bars show (for large 
E and small N) the discrete nature of the data, while vertical error bars show the standard 
error. The solid curve s hows a normal dis tribution with mean 1 and standard deviation 0.3, 
as suggested previously (Singh et al.||2014l; the dashed curve shows the distribution found by 


Boudaoud et al. (20071 from a mean-held aggregation model with a single parameter fitted to 
match experiments. 


of evaporation modifies the dynamics of aggregation in three important ways: (i) any 
clusters that form shrink down more quickly than would be the case without evaporation, 

(ii) clusters are more prone to split up, or bifurcate, when evaporation is included, and 

(iii) the final state generally consists of smaller clusters. The first of these observations 
recalls what was seen in the two-body problem where contact occurs in finite time with 
individual blocks moving at a constant speed, see (3.6). We shall see that the second and 
third observations are intimately related. 

To better understand the third observation, numerical simulations with a range of val¬ 
ues of K and E were performed and the maximum number of blocks in any cluster noted, 
see figure]^ 6). This plot confirms the trend observed in (iii) above. In particular, we note 
that for sufficiently large evaporation rates E, no clusters are formed (the green points in 
figure . Comparing the transition between no clustering and clustering observed with 
many blocks (figure |^) to the transition between separated and stuck states observed in 
the two-body problem (figure [^) we observe a strong similarity; in particular, the line 
KE = 1/2 seems to be important in both transitions. Away from this line, the maximum 
number of elements in a cluster varies as the evaporation rate E varies (at fixed K). To 
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try to understand these results, we begin by considering what happens in the limit of no 
springs. 


4.2. No springs: K = 0 

With no springs, K = 0, the problem simplifies considerably. To see this, we note that 
(2.1) then gives 


Fj+i - 2Fj + Fj_i - 0. 


(4.2) 

If the system is unforced at its edges, Fq = = 0, then this may be inverted to give 

(4.3) 


^ /if dt 


= 0 

hj 


for j = 1,..., iV. The motion within each gap therefore decouples from that in every other 
gap and we recover the two-block problem with K = 0. We conclude that in 

this case all of the blocks will aggregate with contact occurring at / = /* with /» given 
by d^. 


4.3. Finite strength springs K > 0 

Analytical progress is possible in the limit of no springs because the evolution of the 
liquid within each gap becomes decoupled from that in every other gap. However, with 


non-zero spring stiffness, the force balance (2.31 must hold; this couples the evolution 
within the gaps to one another. 

4.3.1. Analytical results close to contact 

To make analytical progress, we consider the limit of small gap thickness, hj 1, so 
that a cluster of blocks are already close to contact. In this limit additional displacements 
are small and the force arising from the springs is approximately constant, simplifying 
the problem considerably. In particular, when the blocks are close to contact, hj I, 


the force-balance equation (2.3) may be approximated by 


-2K = Fj+i - 2Fj 




(4.4) 


Physically, this equation expresses the fact that in the small-gap limit, the net hydro- 
dynamic force on a block, Tj+i — Fj, exceeds that on its neighbour by 2K since it is 
approximately one unit further from its equilibrium position. 


Given two particular Fj somewhere within the system, (4.4) may readily be inverted 


to give an expression for general Fj. Motivated by the repeated cluster-bifurcation events 
seen in figures |^c,c/), we consider whether a cluster of N spring-block elements (enclosing 
N — 1 liquid-filled gaps) that is close to complete collapse will bifurcate. Labelling the 
liquid gaps from one side of the cluster, we can approximate the hydrodynamic force at 
its edges by zero, i.e. Fg = Fjq = 0, since the cluster will be comparatively well separated 


from any neighbouring elements. Inverting (4.4) subject to these edge conditions, we find 
that 

F,=Kj{N-j) (4.5) 

Thus the hydrodynamic force Fj in a gap is given by a constant value /C(j; N) = Kj{N — 
j) that depends only on the label j of that gap within a particular cluster, and the 
number of blocks N in the cluster. (The detailed gap thicknesses do not appear because 
the displacements of blocks from their equilibrium positions, and the consequent spring 
forces, are nearly constant when the blocks are close to contact.) 


Now, recalling from (2.1) that the hydrodynamic force K is the sum of lubrication 
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and capillary components, we can write (4.5) together with the equation describing mass 
conservation as the system 


^3 ‘^^3 

hi d( 


dt 


-^ + ic{r,N) 

h3 

dt hi 


(4.6) 

(4.7) 


These equations are equivalent to those describing the evolution of a single gap, (3.1)- 


(3.2), in the limit of h <C 1, but with an effective spring stiffness lC{j;N). Hence the 


many-body nature of the cluster problem enters only through the effective spring stiffness 
/C(j;iV), which varies through the cluster from relatively soft springs at the edge to 
relatively stiff springs in the centre. Within the cluster the presence of nearby elements 
increases the spring force acting to open an individual gap. This is reminiscent of the 
‘collaborative stiffening’ that occurs when elastic plates are stuck together by surface 
tension: clusters interact just as individual elements do, albeit with a larger effective 


stiffness that comes from having to bend many plates (Bico et al. 2004). Collaborative 


stiffening here corresponds to a multi-block cluster whose springs have total stiffness 
^ NK and typical displacement proportional to N. Hence the typical force in the cluster 
^ N^K. In the two-body problem the displacement close to contact is unity and so, to 
achieve the same force, the effective stiffness satisfies K. ~ N'^K, as seen in (4.5) and 


(4.6). 


In (3.1 we showed for the two-body problem that ii K < Kcrit{E) the gap thickness 


will shrink to zero, and the blocks stick; otherwise the blocks will separate and return 


to their original positions. From the similarity of (4.6)-(4.7) to (3.1)-(3.2) at small gap 
thicknesses, it seems reasonable to expect that, for a cluster containing N blocks to 
collapse completely to a final stuck cluster, the effective stiffness JC{j\N) must be less 
than Kcrit{E) for each of the gaps within that cluster, i.e. 


max{/C(j;iV)} < Kcrit{E). 


(4.8) 


(We assume that any differences between the systems at large gap thicknesses have no 
significant effect on Kcnt{E).) 


Noting that maxj{IC{j; N)} = KN'^/A, we can then restate the criterion in (4.8) as 
collapse of an Wblock cluster (for N ^ 2) can only occur if 


N < IVmax = 




K 


(4.9) 


Thus if iVmax < 2, which occurs at sufhciently large spring stiffnesses, all blocks end 
up separate (in a trivial = 1 cluster). More generally, (4.9) is a prediction for the 


maximum number of blocks in a cluster as a function of the evaporation rate E and 
spring stiffness K. The only assumption used in its derivation was that whether a cluster 
collapses, and the blocks stick, or splits up should be determined by the behaviour of the 
system with small gap thicknesses, hj ^ 1. 

4.3.2. Reassessment of numerical results 

Figure replots data for the maximum number of blocks in any cluster from figure 
i&), at seven different evaporation rates E and with a variety of spring stiffnesses K 
(i.e. the values of A^max along seven different horizontal lines in hgure[^). In figure]^ this 


data has been rescaled based on the analytical prediction (4.9): figureWo) highlights how 
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Figure 6. Rescaling of the data in figure|^6) to test the theoretical prediction (4.9| for the 
maximum number of blocks in a cluster, (a) The maximum number of blocks in any cluster, 


A^max, as a function of Kcv\t{E)/K (points); equation | |4.9[ ) suggests that, provided A^max > 1, 
this should follow the behaviour y = (solid line), (b) An alternative rescaling showing the 

dependence on the evaporation rate, E, directly for data with A^max ^ 2. Here, the solid curve 
sh ows 2KcTit(Ey^^, the dashed horizontal line shows A^maxA^'^^ = 2^/^7r/9 « 3.95 (predicted 
by Singh et al 12014 in the limit E — 0) and the dotted line shows = (2/i?)'^^^, valid 

for A > 0.15. In both panels, colours and symbols correspond to results obtained at different 
evaporation rates: E = 10“^ (red squares), E ~ 10“^ (green triangles), i? ~ 0.1 (blue circles), 
E ~ 1 (black inverted triangles), E = 5 (magenta ‘x’), i? = 10 (cyan ‘+’) and E = 10^ (yellow 
side triangles). 


the raw values of Af^ax vary with Kcnt{E)/K while figure[^6) highlights the dependence 
on E. 

We note first that both rescalings in figure |^ead to reasonable collapses of the data 
from a two-dimensional parameter space (figurel^) onto a single master curve. Since both 
of these rescalings follow naturally from (4.9), we conclude that the reduction from the 
many-body problem to the two-body problem yields useful insight. However, we also note 
that for small evaporation rates, i? <C 1, the maximum cluster size A^max lies somewhat 
below that predicted on the basis of (4.9). Instead, it seems that for Al <C 1 the results 
are closer to the prediction from the no-evaporation problem, E = 0, ( Singh et al.||2014 ) 
where 

(4.10) 


fj — 


3.95. 


Equation (4.10) was obtained from an analysis of front propagation in an unstable 
medium, which was found to leave behind clusters of this well-defined size. Though larger 
clusters would, in principle, be able to collapse and not bifurcate, they are apparently 
not formed by the initial front-propagation mechanism. For all values of E, we find that 
clusters do not reach the maximum size predicted by the linear instability of a nearly 
uniform state without evaporation, which in this notation reads A^maxAT^^^ = 2?! (Singh 

eTZIIMIl ). 


5. Discussion 

We have studied the effect of evaporation on a simple model of elastoeapillary aggre¬ 
gation, which revealed two important differences from the same model in the absence of 
evaporation. Firstly, evaporation acts to bring the elements together more quickly and 
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leads to contact within a finite time. Secondly, the maximum number of elements in a 
final cluster decreases with increasing evaporation rate. 

The first of these observations is important for applications since it suggests that 
evaporation can drive objects together more quickly than would be the case without 


evaporation, thereby increasing the risk of, for example, MEMS stiction (Tanaka et al. 


1993). However, the practical consequences may be (partially) mitigated by the second 


observation: if high evaporation rates cause fewer elements to come into contact, then 
those elements that do come into contact are less likely to be deformed enough to fracture 
(which is one of the main manufacturing concerns, see figure [^). 

Our analysis of the mathematical model allowed us to reduce the problem of a cluster 
that is close to collapsed (i.e. small gap thicknesses) to a two-body problem with an 
effective spring stiffness. If a cluster starts to form then this effective stiffness is smallest 
at the edges and largest at the centre of the cluster; thus, if the cluster breaks up (because 
it is too large) then we expect a void to form in the centre. This is similar to observations 
of microscopic pillars made up of many carbon nanotubes: if the diameter of the initial 
pillar (cluster) is small enough then the entire pillar will collapse during evaporation; if 
not, then voids form in the centre ( de Voider fc Hart]|2013 ). 

We have shown that the number of elements in a cluster is controlled by the dynamics 
of evaporation. As such, the final pattern observed depends on the dynamics of pattern 
formation and not simply which pattern is the global energy minimizer, as is often as¬ 
sumed. However, we have also seen that the effect of evaporation only becomes significant 
once A > 0 . 1 , i.e. when the evaporation dynamics occur on a time scale comparable to, 


or faster than, that of the lubrication-type dynamics itself. From (2.6) the dimensionless 


evaporation rate E = 2Ca,e{xo/w), where Cag = is an evaporative capillary number. 
Written in this way, we see that E is the product of a property of the liquid filling the 
gaps (Cae) and the initial aspect ratio of the liquid gaps {xq/w). For the evaporation of 
pure, still water at room temperature and ambient humidity, e « 60 nm/s ( Li et aZ.|2012 
Routh|[2013 ), fi = 10“^ Pas, 7 = 72 mN/m so that Cag « 10“®. Other liquids evaporate 


more quickly and hence may have a larger value of Cug. For example, HFE-7100 (3M) 
has an evaporation rate e ~ 20 /rm/s (see Lyulin & Kabov 2013 Machrafi et al. 2013 
for example), 7 = 13.6 mN/m and p, = 6.1 x 10“^ Pas so that Gag « 10“^^. Such small 
values of the evaporative capillary number mean that to obtain A > 0.1 would require 
xq/w > 10'^. While this is at least consistent with the requirement xq/w ^ 1 for lubri¬ 
cation theory to be valid, it is unlikely ever to be the case for the micron-scale plates 
and pillars often considered experimentally (see Wei et aTjpOlS for example). However, 
such large aspect ratios are possible for nanotube forests where the typical length of a 
nanotube is on the order of 100 /im but the radius, and hence the typical separation, is 
0(10 nm) ( [Chakrapani et aT]|20(34 ). 

The geometrical enhancement of the evaporation rate suggests that our prediction 
that clusters typically contain fewer elements as the evaporation rate increases may be 
observable in sufficiently dense ‘forests’. Indeed, from their experiments on the formation 
of cells in carbon nanotube forests (see figure [l|)), |Chakrapani et al. ( |2004 ) report that 
“[a faster rate of evaporation] favors crack growth and results in a decrease in the average 
cell width” (i.e. fewer elements form in the clusters between cracks). However, they do 
not present quantitative data to support this comment. 

We hope that the theoretical results presented here may motivate experimental efforts 
to better understand how evaporation rate influences pattern formation in these systems. 
At the same time, there are a number of improvements that should be made to the model 
developed here. For example, experiments are bound to take place in a three-dimensional 


setting where the force from surface tension will require a more detailed calculation (Wei 
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et aL|20f5 1 and the resistive forces from the lubricating flow may be reduced since liquid 


can flow horizontally, as well as vertically. 

Even in situations where the delicate balance between evaporation and lubrication 
dynamics discussed here is not appropriate, other, slower timescales exist with which 
evaporation can more readily compete. For example the viscoelastic relaxation time of the 
elements themselves may be important ( Wei et ar]|2015 |. We expect that the techniques 
developed in this paper may be applied, with suitable modifications, to understand the 
interplay between evaporation and dynamics in such systems. 
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Appendix A. Derivation of model 

In this Appendix, we outline the derivation of the expression for the force exerted 
by the j-th liquid gap on each of its neighbouring blocks, (2.1). More details are given 
by Singh et al. (2014). The well-known lubrication equation, together with the local 
conservation of mass leads to a relationship between the pressure profile within the gap, 
Pj{x), and the rate at which the gap thickness, /ij, is changing with time. In particular, 


dt 12pL dx^ 


dhj 


(Al) 


This equation for pj is to be solved subject to the boundary conditions of no flux through 
the base, i.e. dpj/dx\x=Q = 0, and that the pressure just beneath the meniscus is equal 
to the capillary pressure, i.e. pj(a;“) = — 27 /hj . We find that 


, dp dhj 

= T3- 


dt 




27 


The force on each of the blocks that bound the j-th gap is then given by 




Pj 


dx = 


4px' 


3 

dt 


2'jXj 


(A 2) 


(A3) 


which is the dimensional version of (2.1). 


Appendix B. Phase-plane analysis to determine Kcnt{E) 

ForO < KE < i the phase plane of (3.3)-(3.4| is topologically equivalent to figure 2(a): 
trajectories generically tend either to the stable node (0,1) (separated blocks) or to the 
stable node (^*,0) (stuck blocks). The dividing trajectory between these outcomes is the 
stable manifold of the saddle point (^g, 0), which is given by ^ ^ £,s + ‘^h^s/{^s — ^) + 0{h^) 
as h —>■ 0. Starting with this limit, we can integrate back along the stable manifold to 
h = 1, which defines a point 1). Then if AT ^ Alcrit an initial condition , 1) 

leads to a stuck state. 

As KE Z' i the points (Cs,0) and (^*,0) merge into a saddle-node at (2,0) and 
ATcrit = 3.30 -|- 0(1 — 2KEY^‘^. For KE > ^ all trajectories lead to the separated state. 

At KE = 0 the (nontrivial) nullclines of (3.3)-(3.4) coincide as a line of fixed points 
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h + ^ = 1, corresponding to a mechanical equilibrium K{\ — h) = x/h between the 
spring and capillary forces in the absence of evaporation. Away from this equilibrium, 
the trajectories satisfy = const., corresponding to mass conservation hx = const.. 

To consider the useful limit of slow but non-zero evaporation, we lei e = KE and 
note that the stable manifold of the saddle point, = 1 -|- e/2 — h{l + 3e/4) -|- 0(e^, h^), 
starts almost parallel to the two nearly coincident nullclines. This motivates setting 
Tj = {h + ^ — \)/e so that ( |3.3[ )-( [3^ become 

^ = eC [27? - ^ 2 ] ^ = 

Since e <C 1, we expect the evolution of ^ to be slow and hence, integrating backwards, 
7? collapses onto a ‘slow’ manifold given by d7?/ds = 0(e), or 


7 ? 


3C- 1’ 


(B2) 


over the range 1/3 < ^ < 1. As ^ approaches 1/3, 77 diverges and we need to rescale 
the equations to understand how the stable manifold leaves the 0{e) neighbourhood of 
h + i = l. 

A suitable rescaling is to set x = 3(3^ — l)/2e^/^, y = = 9(^ + h — l)/e^/^. At 

leading order (B 11 becomes 

and thus ‘^ = 2x -. (B 3) 

dx y 


3 dx 
ds 


= y, 


3 dy 


For the solution to (B3) to match onto (B2), we require 7 /^l/ 2 a:asx—>- 00 . This 


condition specifies a unique solution to (B3), given implicitly by a; = —Ai'(z)/Ai( 2 :), 
where Ai denotes the Airy function and z = x'^ — y. As x ^ — 00 , y ^ x^ — Zq, where 
zq = —2.338 is the first zero of Ai(z). This asymptotic behaviour matches onto the 
trajectory = (4/27)(l —e^/^2;o/3) + 0(e), which continues to 1) . Hence ATcrit 

(27/4)(l -I- €^'’’^ 2 : 0 / 3 ) ~ (27/4) — 18.The limiting value K„it = 27/4 for A = 0 was 
obtained by |Singh et al. ( 2014[ ). 
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